word文檔可自由復(fù)制編輯word文檔可自由復(fù)制編輯太陽影子定位模型摘要趨勢；對桿長呈正比關(guān)系增長；對經(jīng)度呈現(xiàn)先急劇增長到峰值再突變?yōu)?1中提供的數(shù)據(jù)，畫出了天安門廣場上直桿的太陽影子分布曲線圖。20 21? i i

? 歸i 歸i|DA-DA'|,min |S

-S'|。i=1

i=1)桿長為，太陽方向角殘差比為1.8%，影長殘差比為型求解的時間復(fù)雜度較高。為提高計算速度，引入了粒子群算法。分別對附件2和3中的數(shù)據(jù)進行分析，確定出的地點坐標(biāo)分別為)，)和))2為西藏阿里，日期為8/14或4/，附件3為廣西梧州市，日期為12/27或12/14?？梢詮囊曨l中直接得到的影長并不是實際長度，而是其投影長度，這里采用基于Hough變換和透視變換的圖像矯正法，對斜視圖像進行矯正，得出實際影長。然后將得到的數(shù)據(jù)帶入問題二的模型中，給出視頻拍攝地點為)題三的基于粒子群算法的目標(biāo)規(guī)劃模型，求解出視頻拍攝地點為)——內(nèi)蒙包頭市境內(nèi)，拍攝日期為6/11或7/13。建筑物群合理間距問題。關(guān)鍵詞：雙目標(biāo)規(guī)劃粒子群算法Hough變換透視變換一、 問題的重述1.1問題的背景與太陽活動有著密切的聯(lián)系。而古代智慧的先民就利用了這個現(xiàn)象制作了日晷，時間和地理信息，是非常方便可靠的。1.2問題的提出的地點和日期的一種方法。針對上述背景和應(yīng)用需求，提出以下問題：1.建立影子長度變化的數(shù)學(xué)模型，分析影子長度關(guān)于各個參數(shù)的變化規(guī)律，并應(yīng)用你們建立的模型畫出2015年10月22日北京時間9:00-15:00之間天安39度54分26秒,東經(jīng)116度23分29米高的直桿的太陽影子長度的變化曲線。1個可能的地點。3.根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型2和附件3的影子頂點坐標(biāo)數(shù)據(jù)，給出若干個可能的地點與日期。4．如果已有一根直桿在太陽下的影子變化的視頻，并且已通過某種方式估計出直桿的高度為2與日期。二、 問題的分析2.1問題一的分析過程。2.2問題二的分析值之和達到最小時，所得經(jīng)緯度即為桿子的地點坐標(biāo)。2.3問題三的分析與問題二不同的是，該問中日期是個未知量。首先考慮沿用上一問的模型。365理的適應(yīng)度函數(shù)，期望提高計算效率。2.4問題四的分析得出實際影長。然后再分別利用問題二和問題三的模型進行求解。三、 問題的假設(shè)1.2.地球是一個完美的球形，不考慮海拔、地球扁率的影響。3.無光線衍射造成的影子減淡現(xiàn)象。5.在小尺度考慮直桿投影問題時，地表是絕對水平的。6.不考慮地球公轉(zhuǎn)的影響。7.題目所給的數(shù)據(jù)是真實的，可靠的。四、 符號說明符號說明q太陽高度角L水平地面上直桿長度S水平地面上直桿影子長度h當(dāng)?shù)氐胤綍r時角d當(dāng)?shù)靥柍嗑暜?dāng)?shù)鼐暥戎礸當(dāng)?shù)亟?jīng)度值t當(dāng)?shù)氐胤綍rt0北京時間n當(dāng)日日期序號P當(dāng)?shù)亟?jīng)緯度坐標(biāo)A當(dāng)?shù)靥柛叨冉荰時間（包含時刻與日期）b離散化后的方位朝向ei每組數(shù)據(jù)影長與x坐標(biāo)軸的夾角xi，yi題目附件提供的第i組直桿坐標(biāo)d日期注：其他符號將在下文中給出具體說明。五、 模型的建立與求解5.1問題一5.1.1模型的建立一些地理學(xué)定義，以及重新定義一些本模型需要用到的參數(shù)。上，直桿長度與影長的比值即為太陽高度角的的正切值： tanq=L S圖1太陽方位角示意圖其中L為桿長，S為影長。通過查閱，得知太陽高度角的計算公式為： sinq=coshcosdcosf+sindsinf (2)其中q為太陽高度角，h為地方時時角，d為當(dāng)時的太陽赤緯，f為當(dāng)?shù)鼐暥?。以一個地方太陽升到最高的地方的時間為正午1212時之間等分為24個小時，所成的時間系統(tǒng)，稱為地方時。地球上每一個地點都有角時，要先將北京時間換算為當(dāng)?shù)氐胤綍rt：g-120?t=t0+

()其中t0為北京時間，g為當(dāng)?shù)亟?jīng)度。根據(jù)某地地方時，可以換算出當(dāng)?shù)氐牡胤綍r時角。地方時時角h即為當(dāng)?shù)嘏c子午線之間相差的角度：h=15按(t-12) (4)緯的最高度數(shù)為'，夏至?xí)r太陽的赤緯為，冬至?xí)r太陽的赤緯為-23?26'。春分和秋分時太陽的赤緯為0?。d可用下面這個公式來計算：d=23.45?sin?(284+n)? (5)? 365 ?? ?其中n為當(dāng)日日期序號，1月1日時n=1，以此類推得10月22日n=聯(lián)立式子(2)-(5)得到方程組：? Ltanq=S?S??sinq=coshcosdcosf+sindsinf?t=t0+?

g-120?15?

(6)?h=15按(t-12)? ? + ?=23.45?sin?(284 n)?求解上述方程組，得：

? 365 ?S=L/tan(arcsin(cos(15按(t-12))cos(d)cosf+sin(d)sinf)) (7)可見，影子長度的變化與當(dāng)?shù)氐乩砦恢肞(j,f)，直桿長度L，時間T(t,d)這三個參數(shù)有關(guān)。5.1.2模型的求解首先計算題目所給條件下的h與f(2)9:00-15:00制成下表以供參考：表1影子長度分布北京時間太陽高度角/度影長/m9:0021.187.7412:0037.883.8512:1437.993.8415:0025.356.33將影長隨時間變化的情況用繪制成圖像：圖2影子長度變化曲線2中已經(jīng)給出；對于地點，將其分為經(jīng)度g和緯度f這兩個參數(shù)進行運算求解，求解后的變化規(guī)律圖如下：圖3影長與日期關(guān)系曲線圖4影長與桿長關(guān)系曲線圖5影長與經(jīng)度關(guān)系曲線圖6影長與緯度關(guān)系曲線5.1.3結(jié)果分析從圖3最小點出現(xiàn)在12:14附近時，會產(chǎn)生時差，顯然是符合常理的。從圖4tanq=LS了模型的準(zhǔn)確性。從圖5地平線以下，就不存在影子，影子長度也為0。從圖65緯度的增加太陽高度角減小，從而導(dǎo)致了影長的增加。5.2問題二——基于網(wǎng)格逼近算法的雙目標(biāo)規(guī)劃模型5.2.1數(shù)據(jù)預(yù)處理x軸的方位朝向角為b(0??b<360?)軸的方位朝向角為0?。21取部分呈現(xiàn)在下表中：表3影子方位角數(shù)據(jù)表北京時間x坐標(biāo)/米y坐標(biāo)/米方位角/度14:421.03560.4973b+25.63?14:571.20870.5255b+23.49?15:121.39550.5541b+21.65?15:271.60330.5833b+20.02?15:421.82770.6135b+18.55?5.2.2模型的建立定義為從北方沿著地平線順時針量度的角：cosA=sind?cosf-cosh?cosd?sinfcosa

(8)上述公式可以用來計算近似的太陽方位角，不過因為公式是使用余弦函數(shù)，180180(下午)，方位角應(yīng)該大于180度，即要取補角的值，故作如下修正：?arccos?sind?cosf-cosh?cosd?sinf?,

h<0? ? cosa ??A=? ? ???360?-arccos?sind?cosf-cosh?cosd?sinf?,h?0

(9)? ? cosa ?? ? ?由上述推導(dǎo)可得，方位角A與坐標(biāo)系方位朝向b的關(guān)系為：=b+ei,

i=,21 (10)其中ei為每組數(shù)據(jù)影長與x坐標(biāo)軸的夾角：ie=arctanyi,i

i=2, ,21 將上一問題的某些參數(shù)關(guān)系與式-列成方程組：??ei?

=arctanyi,

i=2, ,21?=b+ei,?

i=2, ,21? + ?? ?=23.45? ?(284n)365 ?

12)? ? ??h=15按(t-12)?=cos(15按(t-12))cos(d)cosf+sin(d)sinf?t=t0t=t0+? 15?將上述求太陽高度角的方程化為如下關(guān)系式：cos(A)=

(sin(d)?cosf-cos?按(t+15 15 0

g-120?15?

-12)??cos(d)?sinf)

? ? g-120? ??cos(arcsin(cos?15按?t0+ -12??cos(d)cosf+sin(d)sinf))? ? ??下面對各個參量進行影響因素的分析：g-120?t=t0+

，h=15按(t-12)?(284+n)?當(dāng)?shù)亟?jīng)度g的變化以及時間t0有關(guān)；對于赤緯角，由d=23.45?sin? ?? 365 ?為赤緯角dsinq=coshcosdcosf+sindsinf式子可關(guān)，我們可以將其理解為地點的變化；對于太陽方位角，由式子(8)可知，太陽大小有關(guān)，同樣可以將其理解為地點參數(shù)的變化。綜上，我們對參量的影響因素進行總結(jié)，可以得到下表：表4影響因素分析參量影響因素時角(h)經(jīng)度g，時間t0赤緯角)常量太陽高度角)地點坐標(biāo)(g,f)，時間t0太陽方位角(A)地點坐標(biāo)(g,f)，時間t0可以發(fā)現(xiàn)，在式中，方程右端的未知量只有地點坐標(biāo)(g,f)和北京時間t0，故可以將其簡化為如下形式：F,f,t0)=(14)式中，映射F表示方程右端關(guān)于地點坐標(biāo)(g,f)這兩個自變量的函數(shù)。word文檔可自由復(fù)制編輯21個不同時刻對應(yīng)的太陽方位角(i=2,...,21個不同的方程，將其列成方程組如下所示：?F,f,t0)=cos()??F,f,t0)=cos()??,f,t0)=cos()由式(9)(g,f)計算得出此處不同時刻的太陽方位角(7)' 'i據(jù)地點坐標(biāo)(g,f)算出此處不同時刻下的影子長度S'，對其進行歸一化處理，可i以得到：

'Si-'

1211'?Si'S='S=歸i S'

21i=1S'

16)imax

imini對此我們可采用遍歷算法對不同經(jīng)緯度(g,f)下的太陽方位角的插值DA'以i及歸一化后的影子長度S歸i優(yōu)解。

進行求解，從而與實際值作差進行比較，進而得出最故此，我們建立多目標(biāo)規(guī)劃模型：?20?? ? i i?目標(biāo)函數(shù)：??|-|i=121

17)?min |S

-S'|? ? i=1?-180??g?180??約束條件：?-90??f?90????0??q?90??

(18)式中，表示實際的太陽方位角，S歸i表示歸一化后的實際影長。5.2.1模型的求解0.01，對其進行編程求解：表5定點坐標(biāo)數(shù)據(jù)經(jīng)度緯度桿高19.12?表6誤差分析表太陽方向角殘差影長殘差太陽方向角殘差比影長殘差比0.0260.014word文檔可自由復(fù)制編輯即桿子所處地點坐標(biāo)為)，在海南東方市境內(nèi)。5.2.1結(jié)果分析)的影子長度殘差比為，太陽方向角殘差比為1.8%明得到的經(jīng)緯度坐標(biāo))是較為精確的。5.3問題三5.3.1模型的建立變的情況下，要求模型輸出桿子所處的地點和日期。由表3可以看出，日期影響的參量為赤緯角，其他參量并不受日期的影響，太陽赤緯公式為：d=23.45?sin?(284+n)? (19)? 365 ?? ?式中，n為當(dāng)日日期序號，1月1日時n=1，以此類推，日期每增加一天，對n進行加1即可。n及經(jīng)緯度坐標(biāo)(g,f)進行遍歷求解，建立與上一問相似的多目標(biāo)規(guī)劃模型，從而求出最優(yōu)解：?20?? ? i i?目標(biāo)函數(shù)：??|-|i=121

(0)?min |S

-S'|? ? i=1?-180??g?180?約束條件：??-90??f?約束條件：???0??q?90??>0

(21)式中，表示實際的太陽方位角，S歸i表示歸一化后的實際影長。5.3.2模型的求解編程求解：表7桿子所處地點和日期及誤差分析數(shù)據(jù)表附件二附件三經(jīng)度經(jīng)度緯度緯度桿高桿高太陽方向角殘差比太陽方向角殘差比影長殘差比影長殘差比日期8/14或4/日期/27或12/145.3.3模型的改進——基于粒子群算法的目標(biāo)規(guī)劃模型對目標(biāo)函數(shù)進行求解，從而降低模型求解的時間復(fù)雜度。假定有一個D維的目標(biāo)搜索空間有n個微粒組成了一個粒子群其中每個微粒都用一個D維的向量描述，將它的空間位置表示為i=(1,i2,,D),i=,,,n這可看做目標(biāo)優(yōu)化問題中的一個解代入適應(yīng)度函數(shù)計算出適應(yīng)度值可以衡量微粒的優(yōu)劣；第i個微粒的飛行速度也是一個D維的向量，記為i=(1,i2,,D)；第i個微粒所經(jīng)歷過的具有最好適應(yīng)值的位置稱為個體歷史最好位置記為i=(1,i2,,D)整個微粒所經(jīng)歷過的最好位置稱為全局歷史最好位置，記為g=(g1,g2,,D)，粒子群的進化方程可描述為(t+1)=(t)+(t)((t)-(t))+(t)(pgj(t)-(t))(t+1)=(t)+(t+1)

(2)其中i表示第i個微粒，j表示微粒的第j個維度，t表示第t代，c1，c2個加速常量通常的取值范圍是(,2)，1的隨機函。

U(0,1)

U(0,1)是兩個相互獨立可以調(diào)節(jié)微粒去往自身周邊的最好位置，c2可以調(diào)節(jié)微粒去往整個粒子群所能找到的最好位置。散的坐標(biāo)值點：mi=(g,f,d)“粒子Zj了能夠?qū)ふ业阶銐蚝玫膬?yōu)化解，而且避免計算量的增大，本模型中設(shè)置了20個“粒子”來進行優(yōu)化解搜索，即令i=。由于本問中直桿坐標(biāo)和日期都是未知的，因此將粒子設(shè)為儲存三個參數(shù)：?Z1=經(jīng)度g???Z2=緯度f (23)??Z3=日期dg=rand(-180,+180)，f=rand(-90,+90)，d=rand“最優(yōu)解”很可能會是局部最優(yōu)解，而非全局最優(yōu)解。度值Q做定義為：? ?1?2?R=? ??Q=? ?V??

(24)=

20 21|DA-DA'|+ |S

-S'|?? i=1

i i ? 歸i 歸ii=1記錄做比較。如果能夠比最好記錄更好就將此值設(shè)為最好記錄。接下來，參照式(22)，對每個粒子的去向進行調(diào)整。每次調(diào)整時，都要對?(0,2)，?(0,2)，?(0,1)，?(0,1)進行重新隨機取值，保證粒子行為的隨機性：既能有向最優(yōu)解靠近的趨勢，也能存在尋找其他更優(yōu)解的可能。時間尋找解。這里令迭代次數(shù)為1000。5.3.4粒子群算法模型的求解表8桿子所處地點和日期及誤差分析數(shù)據(jù)表附件二附件三經(jīng)度經(jīng)度緯度緯度桿高桿高太陽方向角殘差比太陽方向角殘差比影長殘差比影長殘差比日期8/14或4/日期/27或12/145.3.5模型的綜合結(jié)果分析件12到的日期分別關(guān)于夏至日6/22及冬至日12/對稱，這是由于太陽的運動是關(guān)知參量增加時，可以使用粒子群算法進行求解。5.4問題四5.4.1數(shù)據(jù)的采集及預(yù)處理——基于Hough變換及透視變換的圖像矯正過的相關(guān)知識從圖片中還原出影子的真實長度及變換的角度。透視投影的一個主要特征是可以延伸到無窮遠的物體在照片平面中可以獲點方向。而消隱點在水平方向上的延長即為地平線。據(jù)此，可以引出地平線的性質(zhì)：在透視圖里，從地平線上任一點引出的兩條直線都是平行的。到消隱點，消隱點在水平方向的延長線就是地平線。令地平線與直桿的交點為P則P點就是直桿向其正前方的消隱點從直桿底部引出一水平輔助線B,則從P點到直桿底部O的連線L與水平線垂直。將影子端點R與PAP點引出的A平行于LA與B也垂直。因此，?RQO為直角，三角形DRQO為直角三角形。下面利用Hough變換[5]提取平行直線，我們知道，當(dāng)圖像中存在很多條直線若直線的寬度大于1，則舍棄。這個過程可以表示成：

m=ax+by+c,ox+py+1

(5)n=dx+ey+f (5)式中，斜視圖像點的像素坐標(biāo)為(x,y)，正視圖像的像素坐標(biāo)為(m,n)，a,b,c,d,e,f,o,p,q為透視變換參數(shù)。先，利用對其進行圖像提?。ㄌ崛∷俣葹槊?產(chǎn)生的22幅圖片進行灰度處理和高斯濾波處理:灰度處理高斯濾波中陰影也變得更加清晰。5.4.2模型的建立1.在拍攝時間已知的情況下：代入到問題二的多目標(biāo)規(guī)劃模型中，即可得到視頻拍攝地點(g,f)。?20?? ? i i? 目標(biāo)函數(shù)：??|-|i=121

(6)?min |S

-S'|? ? i=1?-180??g?180?? 約束條件：?-90??f?90????0??q?90??

(27)式中，表示實際的太陽方位角，S歸i表示歸一化后的實際影長。2.在拍攝時間未知的情況下：間復(fù)雜度，并且誤差比很小。題三的目標(biāo)規(guī)劃模型中利用粒子群算法進行求解：?20?? ? i i? 目標(biāo)函數(shù)：??|-|i=121

(8)?min |S

-S'|? ? i=1?-180??g?180? 約束條件：??-90??f? 約束條件：???0??q?90??>0

(29)式中，表示實際的太陽方位角，S歸i表示歸一化后的實際影長。5.4.3模型的求解表9桿子所處地點和日期及誤差分析數(shù)據(jù)表已知日期未知日期經(jīng)度經(jīng)度緯度緯度桿高桿高影長殘差比影長殘差比日期——日期6/11或7/135.4.4結(jié)果分析地點坐標(biāo)為))經(jīng)6頻日期為7/13，兩者剛好關(guān)于夏至日6/22對稱，驗證了模型的可靠性。六、 模型的檢驗6.1誤差分析1.標(biāo)準(zhǔn)解與模型所求解的誤差：由于本模型使用了目標(biāo)規(guī)劃方法對問題進發(fā)現(xiàn)其最大值為2.局部最優(yōu)解造成的地理坐標(biāo)誤差：粒子群算法的本身原理會導(dǎo)致其收斂過程中，有時候需要人為干預(yù)來篩去一些的地區(qū)，可能造成冗余解或缺解。3.視頻數(shù)據(jù)預(yù)處理過程中的誤差：利用Hough變換尋找出的平行線并不一定是處在同一水平面上的，因此做透視變換的時候可能會產(chǎn)生誤差。6.2敏感性分析我們是以每2分鐘一幀的速度，從0時刻開始對其提取了22組數(shù)據(jù)。如果將起始時刻進行變換，即可以取1,5,...,39這些時間節(jié)點的視頻數(shù)據(jù)，在已知日期的情況下帶入模型進行求解，從而檢驗視頻數(shù)據(jù)的改變，對模型結(jié)果產(chǎn)生的影響，得到的結(jié)果如下：表10敏感性檢驗對比數(shù)據(jù)原始數(shù)據(jù)檢驗數(shù)據(jù)經(jīng)度經(jīng)度緯度緯度健性。七、 模型推廣可以應(yīng)用到求建筑物群合理間距問題。平行的鉛垂面,其在地面上的投影為兩平行線,當(dāng)所取的軌跡線上的時間差剛好等于要求的必要日照時間數(shù)時,兩幢建筑之間的垂直距離即為合理的間距,這種方法適用于兩幢建筑物等高的情況。假設(shè)樓體1的投影直線方程為=2的投影直線方程為y2=kx2+b，需要獲取的日照時間數(shù)為ts1與投影曲線的交點(x0,y0)? 2+2?型可求解太陽高度角q=arccot??y0L

?h?? ?t=h +12(xt=h s 15?樓間距d=

b =1+k2

y'-kx'。1+k1+k2模型的優(yōu)點：1.中需要對模型代碼進行部分修改就可以進行計算。2.搜索能力，一般不會陷入局部最優(yōu)。3.實用性高。只需要一段描述日影變化過程的視頻以及相應(yīng)的時間就可以求出當(dāng)?shù)氐慕?jīng)緯度坐標(biāo)。模型的缺點：1.只能應(yīng)用在水平地面上。沒有考慮斜坡和坎坷地面上的情況。2.對于拍攝角度不好的視頻，該模型產(chǎn)生的誤差較大。3.對形狀規(guī)格較為復(fù)雜的建筑，本模型不能適用。九、 參考文獻[1][2]袁新生，邵大宏，LINGO和EXCEL在數(shù)學(xué)建模中的應(yīng)用，北京：科學(xué)出版社，2007。[3]謝金星，薛毅，優(yōu)化建模與LINDO/LINGO軟件，北京：清華大學(xué)出版社，[4]程序設(shè)計與應(yīng)用，高等教育出版社，2012。[5]孫榮峰，快速霍夫變換算法，計算機學(xué)報，24(10)：1104-1109，2001。[6]鄭鵬飛，林大鈞，基于影子軌線反求采光效果的技術(shù)研究，華東理工大學(xué)學(xué)報(自然科學(xué)版)，36(3)：458-463，2010。[7]何援軍，透視和透視投影變換，計算機輔助設(shè)計和圖形學(xué)學(xué)報，17(4)：734-739，2005。[8]武琳，基于太陽陰影軌跡的經(jīng)緯度估計技術(shù)研究，天津大學(xué)，2010。[9]葉德意，基于自適應(yīng)變異粒子群算法的分布式電源選址與容量確定，電網(wǎng)技術(shù)，35(6)：155-157。十、 附錄10.1數(shù)據(jù)表格表11太陽高度角與影長數(shù)據(jù)時間太陽高度角影長921.187959167.7392853859.122.085392227.3935161449.222.966792967.0789613349.323.831452266.7918104359.424.678639376.5288782479.525.507602246.2874837049.626.317568086.0653559259.727.107744185.8605606129.827.877318955.6714418429.928.625463295.4965756091029.351332225.33473245310.130.05406695230.732796985.04599397110.331.38664334.91736631310.432.014721014.79825987810.532.616143074.68805864210.633.190024174.58622312410.733.7354854.49228050310.834.251656914.40581626710.934.737686894.326467111135.192742924.25391487811.135.616019514236.006743424.12812394911.336.364179584.07443159911.436.687636934.02662174911.536.97647433.98453738511.637.230106033.94804464211.737.448007373.91703073811.837.629719483.89140223811.937.774853963.8710836071237.883096753.8560160412.137.954211353.84615654912.237.988041333.8414772912.337.9845123.841965126word文檔可自由復(fù)制編輯12.437.943631143.84762142512.537.8654893.85846207812.637.750257273.8745177512.737.598187323.89583435812.837.409607513.92247378912.937.184919883.9545148611336.924596093.99205454613.136.629172764.03520946313.236.299246474.08411768413.335.935468274435.538537974.19986668813.535.109198314.26711188913.634.648229054.34092552713.734.156441154.42159297613.833.6346714.50944049913.933.083774964.6048405811432.504624074.70821816914.131.898099134.82005799414.231.265086084.94091318614.330.606471845.07141548114.429.923140445.21228735114.529.215969595.36435651114.628.485827685.52857336414.727.73357115.70603210914.826.960041975.8979964414.926.166066196.1059310611525.352451896.331540604表12視頻采集數(shù)據(jù)時間節(jié)點影長12.42402402422.39108108132.35783783842.32483483552.31582582662.27654654772.24654654782.20741741792.176066066102.162462462112.121861862word文檔可自由復(fù)制編輯122.105255255132.09039039142.054144144152.027087087162.003063063171.966996997181.94191.912942943201.900900901211.876876877221.85885885910.1源程序第一問程序：clc;clear;shadow1=zeros(1,601);%初始化不同時間段影長sthet=zeros(1,601);%初始化不同時間段太陽高度角的sin值thet=zeros(1,601);%初始化不同時間段太陽高度角n=295;%1月1日到10月22日間的天數(shù)time=[9:0.01:15];%從9點運算到15點，步長為0.01hh=[-48.60833333:0.15:41.39166667];%9點到15點地方恒星時系統(tǒng)下的時角delta=23.45*sin(2*pi*(284+n)/365);%太陽赤緯height=3;%旗桿高度phi=39.9072222222222;%天安門緯度fori=1:601%從9點運算到15點，共運算601個值sthet(i)=cos(h(i)*pi/180)*cos(delta*pi/180)*cos(phi*pi/180)+sin(delta*pi/180)*sin(phi*pi/180);%計算出太陽高度角的sin值thet(i)=asin(sthet(i))*180/pi;%計算出太陽高度角shadow1(i)=sqrt((height/sthet(i))^2-(height)^2);%計算影長endplot(time,shadow1)%輸出圖形ylabel('影長/m');xlabel('時間/時');title('影子長度與時間的關(guān)系');gridon;第二問程序clcclearall%附件一數(shù)據(jù)x=[1.03651.06991.10381.13831.17321.20871.2448...1.28151.31891.35681.39551.43491.47511.516...1.55771.60031.64381.68821.73371.78011.8277];y=[0.49730.50290.50850.51420.51980.52550.5311...0.53680.54260.54830.55410.55980.56570.5715...0.57740.58330.58920.59520.60130.60740.6135];length=sqrt(x.^2+y.^2);sOLen=sum(length);Length=length;length=(length(2:21)-length(1:20))./(length(2)-length(1));angle=acos((x(2:21).*x(1:20)+y(2:21).*y(1:20))./(sqrt(x(2:21).^2+y(2:21).^2).*sqrt(x(1:20).^2+y(1:20).^2)));result=[]; %經(jīng)緯度t=14.7:0.05:15.7; %北京時間min=0.05;kkk=0.1;lon=80:kkk:120;beginl=10*ones(1,size(lon,2));endl=40*ones(1,size(lon,2));n=108;%2015年4月18日delta=23.45*sin(2*pi*(284+n)/365); %赤緯角fori=1:size(lon,2)ts=t-(120-lon(1,i))/15;%當(dāng)?shù)氐胤綍romega=15*(ts-12); %時角beginl(1,i);endl(1,i);forlat=beginl(1,i):kkk:endl(1,i)phi=lat; %緯度h=asin(sin(phi*pi/180)*sin(delta*pi/180)+cos(phi*pi/180)*cos(delta*pi/180)*cos(omega*pi/180));%太陽高度角 a=2*pi-acos((sin(delta*pi/180)-sin(h)*sin(phi*pi/180))./(cos(h)*cos(phi*pi/180)));%太陽方位角l=100./tan(h);%影長l=(l(2:21)-l(1:20))./(l(2)-l(1));a=a(2:21)-a(1:20);% if(sum(abs(l-length))<min)% ganH=mean(Length.*tan(h));% result=[result;lon(1,i)latganHsum(abs(l-length)) n];if(sum(abs(l-length))<min)%尋找最小點min=sum(abs(l-length)); %更新min為兩者更大的ganH=mean(Length.*tan(h));result=[lon(1,i)latganHabs(l-length)abs(a-angle)n];%最小點保存在result內(nèi)endclearh;endn

endiplot(result(:,1),result(:,2),'r*')function[ParSwarm,OptSwarm]=BaseStepPso(ParSwarm,OptSwarm,ParticleScope,MaxW,MinW,LoopCount,CurCount)ifnargin~=7 %輸入容錯error('error')endifnargout~=2 %輸出容錯error('error')end%線形遞減策略w=MaxW-CurCount*((MaxW-MinW)/LoopCount);%得到粒子群群體大小以及一個粒子維數(shù)的信息[ParRow,ParCol]=size(ParSwarm);%得到粒子的維數(shù)ParCol=(ParCol-1)/2;SubTract1=OptSwarm(1:ParRow,:)-ParSwarm(:,1:ParCol);%求解出歷史最優(yōu)值與當(dāng)前位置的差值c1=2;c2=2;%完成一次粒子位置速度最優(yōu)值的更新迭代forrow=1:ParRowSubTract2=OptSwarm(ParRow+1,:)-ParSwarm(row,1:ParCol);%計算出全局最優(yōu)值與當(dāng)前該粒子位置的差值%速度更新公式TempV=w.*ParSwarm(row,ParCol+1:2*ParCol)+c1*unifrnd(0,1).*SubTract1(row,:)+c2*unifrnd(0,1).*SubTract2;%限制速度的代碼forh=1:ParColifTempV(:,h)>ParticleScope(h,2)TempV(:,h)=ParticleScope(h,2);endifTempV(:,h)<-ParticleScope(h,2)TempV(:,h)=-ParticleScope(h,2)+1e-10;%加1e-10防止適應(yīng)度函數(shù)被零除end

end%更新該粒子速度值ParSwarm(row,ParCol+1:2*ParCol)=TempV;a=0.729;%約束因子TempPos=ParSwarm(row,1:ParCol)+a*TempV;%限制位置范圍的代碼forh=1:ParColifTempPos(:,h)>ParticleScope(h,2)TempPos(:,h)=ParticleScope(h,2);被零除

endifTempPos(:,h)<=ParticleScope(h,1)TempPos(:,h)=ParticleScope(h,1)+1e-10;%加1e-10防止適應(yīng)度函數(shù)end

end%更新該粒子位置值ParSwarm(row,1:ParCol)=TempPos;%計算每個粒子的新的適應(yīng)度值ParSwarm(row,2*ParCol+1)=AdaptFunc(ParSwarm(row,1:ParCol));ifParSwarm(row,2*ParCol+1)>AdaptFunc(OptSwarm(row,1:ParCol))OptSwarm(row,1:ParCol)=ParSwarm(row,1:ParCol);end

endword文檔可自由復(fù)制編輯%for循環(huán)結(jié)束%尋找適應(yīng)度函數(shù)值最大的解在矩陣中的位置(行數(shù))，進行全局最優(yōu)值的改變[maxValue,row]=max(ParSwarm(:,2*ParCol+1));ifAdaptFunc(ParSwarm(row,1:ParCol))>AdaptFunc(OptSwarm(ParRow+1,:))OptSwarm(ParRow+1,:)=ParSwarm(row,1:ParCol);end圖像處理程序%測試霍夫變換clcclearcloseall%讀取圖像I =imread('circuit.tif');rotI=imrotate(I,80,'crop');%旋轉(zhuǎn)33度，保持原圖片大小fig1=imshow(rotI);%提取邊BW=edge(rotI,'canny');figure,imshow(BW);%霍夫變換[H,theta,rho]=hough(BW);%計算二值圖像的標(biāo)準(zhǔn)霍夫變換，H為霍夫變換矩陣，theta,rho為計算霍夫變換的角度和半徑值figure,imshow(imadjust(mat2gray(H)),[],'XData',theta,'YData',rho,...'InitialMagnification','fit');xlabel('\theta(degrees)'),ylabel('\rho');axison,axisnormal,holdon;colormap(hot)%顯示霍夫變換矩陣中的極值點P=houghpeaks(H,50,'threshold',ceil(0.3*max(H(:))));%從霍夫變換矩陣H中提取5個極值點x=theta(P(:,2));y=rho(P(:,1));plot(x