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

1、太陽影子定位摘要本文對太陽影子定位問題進(jìn)行了深入探討，建立了太陽影子和時(shí)間以及變位參數(shù)之間關(guān)系的數(shù)學(xué)模型，應(yīng)用了matlab進(jìn)行求解。針對問題一：已知時(shí)間與地點(diǎn)，求解影長與時(shí)間的關(guān)系。通過影長與太陽高度角、太陽方位角的關(guān)系，建立參數(shù)模型，運(yùn)用matlab畫出影子長度變化曲線，影子長度與太陽高度角、太陽方位角、時(shí)間之間的關(guān)系如表1。表1 2015年10月22日9:00-15:00天安門廣場直桿（3米）的太陽影子長度當(dāng)?shù)貢r(shí)間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與時(shí)間t以及影子的觀測坐標(biāo)P，經(jīng)度、緯度、桿長h未知。分析太陽高度角、太陽方位角A與參數(shù)h,的關(guān)系，求解理論坐標(biāo)Q。根據(jù)運(yùn)籌學(xué)理論，以各PQ2之和最小為目標(biāo)函數(shù)，以h,的取值范圍為約束條件，建立非線性規(guī)劃模型。應(yīng)用最小二乘參數(shù)估計(jì)法思想，運(yùn)用窮舉法和放縮法，利用matlab求解，得出,的值。針對問題三，已知t與P，在問題二的基礎(chǔ)上增加一個(gè)未知參數(shù)日期序號n，

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

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

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

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

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

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

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

10、值時(shí) (上午)，方位角的角度小于180度，時(shí)角為正值時(shí) (下午)，方位角應(yīng)該大于180度，即要取補(bǔ)角的值。 cosA=sincos-coswcossincos （8） cosA=sin-sinsincoscos （9）5.1.3建立參數(shù)方程定義圖1中y軸方向?yàn)檎狈较颍白禹旤c(diǎn)坐標(biāo)為x，yx，yL1hz y x A圖1 “立竿見影”坐標(biāo)系由影子長度與太陽高度角的關(guān)系 L1=htan （10）可得出參數(shù)方程： x=sinAL1y=cosAL1 t為參數(shù) （11）根據(jù)建立的參數(shù)方程，運(yùn)用matlab擬合出影子長度與時(shí)間的變化曲線。5.2問題二的模型建立假設(shè)太陽影子頂點(diǎn)坐標(biāo)實(shí)際觀測值為PPxi,Py

11、i,太陽影子頂點(diǎn)坐標(biāo)理論值為Q(QxiQyi),假設(shè)P所在的坐標(biāo)系和正北方向的夾角為0，則有 =hsh, A=hah, Qxi=L1sinAQyi=L1cosA (12) L1=htan (13)由上述三個(gè)式子得到： Qxi=htansinA1Qyi=htancosA1 (14) 化解得到： Qxi=-sinwcoshsin (15) Qyi=sincos-coswcossinsinhQyi=sin-sinsinsincosh （16）注：當(dāng)時(shí)角為負(fù)值時(shí) (上午)，方位角的角度小于180度，時(shí)角為正值時(shí) (下午)，方位角應(yīng)該大于180度，即要取補(bǔ)角的值。通過建立理論坐標(biāo)QQxi，Qyi與觀測坐

12、標(biāo)PPxi，Pyi的距離的平方和S最小的非線性規(guī)模型進(jìn)行求解。即： minS=i=121（Pxi-Qxi）2+（Qxi-Qyi）2 (17)s.t.0360-9090 2h3 5.3問題三的模型建立在第二問的基礎(chǔ)上增加一個(gè)未知參數(shù)n ,建立非線性規(guī)劃模型： min S=i=121Qxi-Pxi2+Qyi-Pyi2 （18）s.t.0360-90902h31n3655.4問題四的模型建立 要根據(jù)一根高度為2米的直桿在太陽下的影子變化的視頻，建立確定視頻拍攝地點(diǎn)。就是要利用視頻幀實(shí)現(xiàn)經(jīng)緯度估計(jì)的方法，能夠根據(jù)未經(jīng)過校準(zhǔn)的影子的位置，得出所拍攝圖像的經(jīng)緯度信息5。具體步驟5如下：（1） 獲取視頻中含

13、有至少兩個(gè)影子軌跡的視頻幀；（2） 對每一幀圖像，檢測影子軌跡；（3） 找出滅點(diǎn)，計(jì)算滅點(diǎn)；（4） 由滅點(diǎn)擬合出地平線；（5） 在圖像沿著有垂直關(guān)系的物體畫出兩條直線，根據(jù)這兩條直線分別與地平線的交點(diǎn)坐標(biāo)，這兩個(gè)坐標(biāo)就是兩個(gè)相互垂直的滅點(diǎn)vx,vy的坐標(biāo)。（6） 根據(jù)公式f=v0vx3vy2+u0vx1vy3+v0vx2vy3+u0vx3vy1-vx1vy1-vx2vy2vx3vy3-u02-v02計(jì)算出焦距f,式中，u0,v0為圖像中心點(diǎn)的坐標(biāo)，根據(jù)公式g=10-u001-v0-u0-v0u02+v02+f2計(jì)算出g（7） 在地平線上找到另一對滿足vxTgvy=0限制的點(diǎn)vp,vq；（8）

14、實(shí)現(xiàn)從圖像到經(jīng)過度量糾正的世界坐標(biāo)的轉(zhuǎn)換；（9） 根據(jù)X=Hx還原出經(jīng)過度量糾正的世界坐標(biāo)中的點(diǎn)，其中X是經(jīng)過度量糾正的世界坐標(biāo)中的點(diǎn)，x是圖像坐標(biāo)中的點(diǎn)；（10） 擬合影子軌跡；（11） 利用日晷原理，計(jì)算緯度；（12） 運(yùn)用極值法求出經(jīng)度。（13） 根據(jù)得到的經(jīng)緯度給出可能的視頻拍攝地點(diǎn)。六、模型的求解6.1問題一的模型的求解曲線擬合根據(jù)參數(shù)方程模型，代入數(shù)據(jù)日期序號n=295，當(dāng)?shù)貣|經(jīng)=116.8667，北緯=39.9072，桿高h(yuǎn)=3米。以y軸為正北方，求解直桿的頂點(diǎn)坐標(biāo)（x，y），得到2015年10月22日北京時(shí)間9:00-15:00之間天安門廣場3米高的直桿的太陽影子的運(yùn)動(dòng)軌跡如圖

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

16、939.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ù)與日常生活中實(shí)際現(xiàn)象相符合，可見此模型是合理的。6.2問題二的模型求解最小二乘法用最小二乘參數(shù)估計(jì)法確定參數(shù)。最小二乘參數(shù)估計(jì)法基本思想：根據(jù)Sh,t的關(guān)系表達(dá)式求出太陽影子的理

17、論坐標(biāo)，通過計(jì)算理論坐標(biāo)與觀測坐標(biāo)的距離，建立距離平方和最小模型。通過對,h進(jìn)行等間距的窮舉最終求得理論坐標(biāo)QQxi，Qyi與觀測坐標(biāo)PPxi，Pyi的距離的平方和S，當(dāng)S取得最小值時(shí)，此時(shí)，h的值即為所求的最佳值。即求解如下最小二乘擬合模型。算法描述：輸入：時(shí)間t、日期序號n輸出：經(jīng)度、緯度、桿高h(yuǎn)的值Step1: 建立以y軸為正北方向的空間直角坐標(biāo)系，以地面為水平面，直桿垂直于地面，建立理論坐標(biāo)QQxi，Qyi的參數(shù)方程，根據(jù)太陽高度角、太陽方位角與影長的關(guān)系，求出理論值QQxi，Qyi關(guān)于，h的關(guān)系式。Step2:時(shí)間t1，t2，t3t21，日期序號n=108，等間距窮舉（附錄5）經(jīng)度0

18、,360，緯度-90,90，桿高h(yuǎn)2，3，根據(jù)理論值QQxi，Qyi關(guān)于，h的關(guān)系式，求得理論坐標(biāo)QQxi，Qyi。Step3:根據(jù)求出的理論坐標(biāo)QQxi，Qyi與附錄1中的觀測坐標(biāo) PiPxi，Pyi （i0,21）,求出距離平方和S=i=121（Pxi-Qxi）2+（Qxi-Qyi）2令S0初始值為，SS0,令S0=S，循環(huán)找出S0的最小值，最小值S=S0。當(dāng)S取得最小值時(shí)，求得，h的值。算法結(jié)束運(yùn)用窮舉法與放縮法，用matlab求解得出最優(yōu)解，如表2。取值范圍0：30：360-90：30：902：0.1：3Step_1所得值120032.9145取值范圍90：5：120-30:5:302

19、：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，h=2.8m。6.3問題三的模型求解在問題二的基礎(chǔ)上，增加一個(gè)未知參數(shù)n，求解如下非線性規(guī)劃模型：minSh,n=i=121（Pxi-Qxi）2+（Qxi-Qyi）2在matlab中根據(jù)n的取值范圍1，365作為限制條件，引入一

20、個(gè)for循環(huán)，同理運(yùn)用窮舉法和放縮法進(jìn)行求解。七、模型的評價(jià)7.1模型的優(yōu)點(diǎn)（1）通過建立非線性規(guī)劃模型,在利用未知參數(shù)取值范圍的條件下，運(yùn)用matlab編程得出了影子頂點(diǎn)坐標(biāo)理論值與實(shí)際值的距離平方之和，最小二乘法求解最小值，該方法結(jié)果合理可信。（2）運(yùn)用matlab繪制出圖形，用數(shù)形結(jié)合的方法來進(jìn)行分析，使模型思路更加清晰，更有說服力。（3）建立參數(shù)方程，可直接求解與逆向求解，模型具有普遍性。（4）通過對參數(shù)范圍的縮小，減小了大量的運(yùn)算，同時(shí)提高了精度，使目標(biāo)函數(shù)的最小值達(dá)到一個(gè)很小的水平，可以認(rèn)為實(shí)際值和理論值幾乎相等S00.0007.7.2模型的缺點(diǎn)（1）由于未考慮實(shí)際觀測值P所在的坐

21、標(biāo)系和正北方向存在的夾角，導(dǎo)致問題二的求解出的值存在一定的誤差。（2）由于運(yùn)用窮舉法，未知參數(shù)過多，等間距窮舉參數(shù)有限，運(yùn)算量過大，存在誤差。（3）模型假設(shè)過于理想化，和實(shí)際情況存在一定的差異。八、模型推廣本文所建的模型可以作為太陽影子定位方面的參考，具有借鑒意義；可以通過該地經(jīng)緯度、時(shí)間、日期、桿高擬合影長隨時(shí)間的變化曲線，可通過該地的影長頂點(diǎn)坐標(biāo)與時(shí)間求解所處位置等?？梢詰?yīng)用在野外求生時(shí)地理位置的預(yù)測、某地的經(jīng)緯度預(yù)測等九、參考文獻(xiàn)1 方榮生，太陽能應(yīng)用技術(shù).北京：中國農(nóng)業(yè)機(jī)械出版社，1985.2 邱國全，太陽時(shí)的計(jì)算.科普資訊，1998. 3 Wikipedia，https:/en.m.

22、/wiki/Solar_zenith_angle，2015,9,144 Wikipedia，/wiki/Solar_azimuth_angle，2015,9,14.5操曉春，曲彥齡，孫濟(jì)洲，武琳，郭曉杰，張煒.基于視頻中太陽影子軌跡的經(jīng)緯度估計(jì)方法.天津：中華人民共和國國家知識產(chǎn)權(quán)局，2009.6地球在線，7韓中庚等.數(shù)學(xué)建模使用教程.北京：高等教育出版社，2012.十、附錄附錄1%problem_1程序如下：clcclear allclose allbjtime = 9:1:15;jingdu = 116+23/60+2

23、9/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:12;w = 15.*(ts1-12)alpha = asind(sind(fai).*sind(sigma)+cos

24、d(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)./(cosd(alpha).*cosd(fai)hold onL1 = 3./tand(alpha)x1 = -sind(A2).*3./tand(alpha);x2 = sind(A2).*3./tand(alpha);y1 = cosd(A1).*3./tan

25、d(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程序如下：clcclear 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*

26、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; w = 15.*(ts1-12)alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w)L=3./tand(alpha)plot(bjtime,L)xlabel(北京時(shí)間)ylabel(影子長度)grid on附錄3時(shí)間、太陽方位角與影子長度的關(guān)系圖：將三維圖旋轉(zhuǎn)一定角度后

27、可得到時(shí)間、太陽方位角與影子長度的關(guān)系圖：太陽方位角與影子長度的關(guān)系圖：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 = 39+54/60+26/3600sigma = 23.45 .*

28、sin(2.*pi*(284+n)/365)ts = 9:1/6:15;w = 15.*(ts-12).*pi./180alpha = asind(sind(fai).*sind(sigma)+cosd(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).*co

29、sd(w)L1 = 3./tand(alpha)plot3(ts,A1,L1,r*-)xlabel(時(shí)刻)ylabel(方位角/度)zlabel(影子高度/米)附錄4時(shí)間、太陽高度角、太陽方位角以及相應(yīng)影子長度的數(shù)據(jù)如下表：時(shí)間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.959

30、396.510756.121615.782065.484465.222884.99264.794.6124.4563時(shí)間10: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.848

31、23.84113.8482時(shí)間12:2012:3012:4012:5013:0013:1013:2013:3013:4013: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時(shí)間14:0014:1014:2014:3014:4014:5015:

32、00高度角31.00129.872928.678627.422426.10824.739123.32A139.619239.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 =

33、 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*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

34、.5368 0.5426 0.5483 0.5541 0.5598 0.5657 0.5715 0.5774 0.5833 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)

35、.*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 smin min = s ; d = jingdu fai h; end end end endmin,d%step_2for jingdu = 90:5:150 ts = bjtim

36、e - 4.*(120-jingdu)+ET; w = 15./(ts-12); for fai= -30:5:30 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 smin min = s ; d =