2015年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽國(guó)家一等獎(jiǎng)?wù)撐?A題-太陽(yáng)影子定位模型的分析

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

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

-S'|。i=1

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

()其中t0為北京時(shí)間，g為當(dāng)?shù)亟?jīng)度。根據(jù)某地地方時(shí)，可以換算出當(dāng)?shù)氐牡胤綍r(shí)時(shí)角。地方時(shí)時(shí)角h即為當(dāng)?shù)嘏c子午線(xiàn)之間相差的角度：h=15按(t-12) (4)緯的最高度數(shù)為'，夏至?xí)r太陽(yáng)的赤緯為，冬至?xí)r太陽(yáng)的赤緯為-23?26'。春分和秋分時(shí)太陽(yáng)的赤緯為0?。d可用下面這個(gè)公式來(lái)計(jì)算：d=23.45?sin?(284+n)? (5)? 365 ?? ?其中n為當(dāng)日日期序號(hào)，1月1日時(shí)n=1，以此類(lèi)推得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)可見(jiàn)，影子長(zhǎng)度的變化與當(dāng)?shù)氐乩砦恢肞(j,f)，直桿長(zhǎng)度L，時(shí)間T(t,d)這三個(gè)參數(shù)有關(guān)。5.1.2模型的求解首先計(jì)算題目所給條件下的h與f(2)9:00-15:00制成下表以供參考：表1影子長(zhǎng)度分布北京時(shí)間太陽(yáng)高度角/度影長(zhǎng)/m9:0021.187.7412:0037.883.8512:1437.993.8415:0025.356.33將影長(zhǎng)隨時(shí)間變化的情況用繪制成圖像：圖2影子長(zhǎng)度變化曲線(xiàn)2中已經(jīng)給出；對(duì)于地點(diǎn)，將其分為經(jīng)度g和緯度f(wàn)這兩個(gè)參數(shù)進(jìn)行運(yùn)算求解，求解后的變化規(guī)律圖如下：圖3影長(zhǎng)與日期關(guān)系曲線(xiàn)圖4影長(zhǎng)與桿長(zhǎng)關(guān)系曲線(xiàn)圖5影長(zhǎng)與經(jīng)度關(guān)系曲線(xiàn)圖6影長(zhǎng)與緯度關(guān)系曲線(xiàn)5.1.3結(jié)果分析從圖3最小點(diǎn)出現(xiàn)在12:14附近時(shí)，會(huì)產(chǎn)生時(shí)差，顯然是符合常理的。從圖4tanq=LS了模型的準(zhǔn)確性。從圖5地平線(xiàn)以下，就不存在影子，影子長(zhǎng)度也為0。從圖65緯度的增加太陽(yáng)高度角減小，從而導(dǎo)致了影長(zhǎng)的增加。5.2問(wèn)題二——基于網(wǎng)格逼近算法的雙目標(biāo)規(guī)劃模型5.2.1數(shù)據(jù)預(yù)處理x軸的方位朝向角為b(0??b<360?)軸的方位朝向角為0?。21取部分呈現(xiàn)在下表中：表3影子方位角數(shù)據(jù)表北京時(shí)間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模型的建立定義為從北方沿著地平線(xiàn)順時(shí)針量度的角：cosA=sind?cosf-cosh?cosd?sinfcosa

(8)上述公式可以用來(lái)計(jì)算近似的太陽(yáng)方位角，不過(guò)因?yàn)楣绞鞘褂糜嘞液瘮?shù)，180180(下午)，方位角應(yīng)該大于180度，即要取補(bǔ)角的值，故作如下修正：?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ù)影長(zhǎng)與x坐標(biāo)軸的夾角：ie=arctanyi,i

i=2, ,21 將上一問(wèn)題的某些參數(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?將上述求太陽(yáng)高度角的方程化為如下關(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))? ? ??下面對(duì)各個(gè)參量進(jìn)行影響因素的分析：g-120?t=t0+

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

'Si-'

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

21i=1S'

16)imax

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

進(jìn)行求解，從而與實(shí)際值作差進(jìn)行比較，進(jìn)而得出最故此，我們建立多目標(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)式中，表示實(shí)際的太陽(yáng)方位角，S歸i表示歸一化后的實(shí)際影長(zhǎng)。5.2.1模型的求解0.01，對(duì)其進(jìn)行編程求解：表5定點(diǎn)坐標(biāo)數(shù)據(jù)經(jīng)度緯度桿高19.12?表6誤差分析表太陽(yáng)方向角殘差影長(zhǎng)殘差太陽(yáng)方向角殘差比影長(zhǎng)殘差比0.0260.014word文檔可自由復(fù)制編輯即桿子所處地點(diǎn)坐標(biāo)為)，在海南東方市境內(nèi)。5.2.1結(jié)果分析)的影子長(zhǎng)度殘差比為，太陽(yáng)方向角殘差比為1.8%明得到的經(jīng)緯度坐標(biāo))是較為精確的。5.3問(wèn)題三5.3.1模型的建立變的情況下，要求模型輸出桿子所處的地點(diǎn)和日期。由表3可以看出，日期影響的參量為赤緯角，其他參量并不受日期的影響，太陽(yáng)赤緯公式為：d=23.45?sin?(284+n)? (19)? 365 ?? ?式中，n為當(dāng)日日期序號(hào)，1月1日時(shí)n=1，以此類(lèi)推，日期每增加一天，對(duì)n進(jìn)行加1即可。n及經(jīng)緯度坐標(biāo)(g,f)進(jìn)行遍歷求解，建立與上一問(wèn)相似的多目標(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)式中，表示實(shí)際的太陽(yáng)方位角，S歸i表示歸一化后的實(shí)際影長(zhǎng)。5.3.2模型的求解編程求解：表7桿子所處地點(diǎn)和日期及誤差分析數(shù)據(jù)表附件二附件三經(jīng)度經(jīng)度緯度緯度桿高桿高太陽(yáng)方向角殘差比太陽(yáng)方向角殘差比影長(zhǎng)殘差比影長(zhǎng)殘差比日期8/14或4/日期/27或12/145.3.3模型的改進(jìn)——基于粒子群算法的目標(biāo)規(guī)劃模型對(duì)目標(biāo)函數(shù)進(jìn)行求解，從而降低模型求解的時(shí)間復(fù)雜度。假定有一個(gè)D維的目標(biāo)搜索空間有n個(gè)微粒組成了一個(gè)粒子群其中每個(gè)微粒都用一個(gè)D維的向量描述，將它的空間位置表示為i=(1,i2,,D),i=,,,n這可看做目標(biāo)優(yōu)化問(wèn)題中的一個(gè)解代入適應(yīng)度函數(shù)計(jì)算出適應(yīng)度值可以衡量微粒的優(yōu)劣；第i個(gè)微粒的飛行速度也是一個(gè)D維的向量，記為i=(1,i2,,D)；第i個(gè)微粒所經(jīng)歷過(guò)的具有最好適應(yīng)值的位置稱(chēng)為個(gè)體歷史最好位置記為i=(1,i2,,D)整個(gè)微粒所經(jīng)歷過(guò)的最好位置稱(chēng)為全局歷史最好位置，記為g=(g1,g2,,D)，粒子群的進(jìn)化方程可描述為(t+1)=(t)+(t)((t)-(t))+(t)(pgj(t)-(t))(t+1)=(t)+(t+1)

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

U(0,1)

U(0,1)是兩個(gè)相互獨(dú)立可以調(diào)節(jié)微粒去往自身周邊的最好位置，c2可以調(diào)節(jié)微粒去往整個(gè)粒子群所能找到的最好位置。散的坐標(biāo)值點(diǎn)：mi=(g,f,d)“粒子Zj了能夠?qū)ふ业阶銐蚝玫膬?yōu)化解，而且避免計(jì)算量的增大，本模型中設(shè)置了20個(gè)“粒子”來(lái)進(jìn)行優(yōu)化解搜索，即令i=。由于本問(wèn)中直桿坐標(biāo)和日期都是未知的，因此將粒子設(shè)為儲(chǔ)存三個(gè)參數(shù)：?Z1=經(jīng)度g???Z2=緯度f(wàn) (23)??Z3=日期dg=rand(-180,+180)，f=rand(-90,+90)，d=rand“最優(yōu)解”很可能會(huì)是局部最優(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è)為最好記錄。接下來(lái)，參照式(22)，對(duì)每個(gè)粒子的去向進(jìn)行調(diào)整。每次調(diào)整時(shí)，都要對(duì)?(0,2)，?(0,2)，?(0,1)，?(0,1)進(jìn)行重新隨機(jī)取值，保證粒子行為的隨機(jī)性：既能有向最優(yōu)解靠近的趨勢(shì)，也能存在尋找其他更優(yōu)解的可能。時(shí)間尋找解。這里令迭代次數(shù)為1000。5.3.4粒子群算法模型的求解表8桿子所處地點(diǎn)和日期及誤差分析數(shù)據(jù)表附件二附件三經(jīng)度經(jīng)度緯度緯度桿高桿高太陽(yáng)方向角殘差比太陽(yáng)方向角殘差比影長(zhǎng)殘差比影長(zhǎng)殘差比日期8/14或4/日期/27或12/145.3.5模型的綜合結(jié)果分析件12到的日期分別關(guān)于夏至日6/22及冬至日12/對(duì)稱(chēng)，這是由于太陽(yáng)的運(yùn)動(dòng)是關(guān)知參量增加時(shí)，可以使用粒子群算法進(jìn)行求解。5.4問(wèn)題四5.4.1數(shù)據(jù)的采集及預(yù)處理——基于Hough變換及透視變換的圖像矯正過(guò)的相關(guān)知識(shí)從圖片中還原出影子的真實(shí)長(zhǎng)度及變換的角度。透視投影的一個(gè)主要特征是可以延伸到無(wú)窮遠(yuǎn)的物體在照片平面中可以獲點(diǎn)方向。而消隱點(diǎn)在水平方向上的延長(zhǎng)即為地平線(xiàn)。據(jù)此，可以引出地平線(xiàn)的性質(zhì)：在透視圖里，從地平線(xiàn)上任一點(diǎn)引出的兩條直線(xiàn)都是平行的。到消隱點(diǎn)，消隱點(diǎn)在水平方向的延長(zhǎng)線(xiàn)就是地平線(xiàn)。令地平線(xiàn)與直桿的交點(diǎn)為P則P點(diǎn)就是直桿向其正前方的消隱點(diǎn)從直桿底部引出一水平輔助線(xiàn)B,則從P點(diǎn)到直桿底部O的連線(xiàn)L與水平線(xiàn)垂直。將影子端點(diǎn)R與PAP點(diǎn)引出的A平行于LA與B也垂直。因此，?RQO為直角，三角形DRQO為直角三角形。下面利用Hough變換[5]提取平行直線(xiàn)，我們知道，當(dāng)圖像中存在很多條直線(xiàn)若直線(xiàn)的寬度大于1，則舍棄。這個(gè)過(guò)程可以表示成：

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

(5)n=dx+ey+f (5)式中，斜視圖像點(diǎn)的像素坐標(biāo)為(x,y)，正視圖像的像素坐標(biāo)為(m,n)，a,b,c,d,e,f,o,p,q為透視變換參數(shù)。先，利用對(duì)其進(jìn)行圖像提?。ㄌ崛∷俣葹槊?產(chǎn)生的22幅圖片進(jìn)行灰度處理和高斯濾波處理:灰度處理高斯濾波中陰影也變得更加清晰。5.4.2模型的建立1.在拍攝時(shí)間已知的情況下：代入到問(wèn)題二的多目標(biāo)規(guī)劃模型中，即可得到視頻拍攝地點(diǎn)(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)式中，表示實(shí)際的太陽(yáng)方位角，S歸i表示歸一化后的實(shí)際影長(zhǎng)。2.在拍攝時(shí)間未知的情況下：間復(fù)雜度，并且誤差比很小。題三的目標(biāo)規(guī)劃模型中利用粒子群算法進(jìn)行求解：?20?? ? i i? 目標(biāo)函數(shù)：??|-|i=121

(8)?min |S

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

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

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

b =1+k2

y'-kx'。1+k1+k2模型的優(yōu)點(diǎn)：1.中需要對(duì)模型代碼進(jìn)行部分修改就可以進(jìn)行計(jì)算。2.搜索能力，一般不會(huì)陷入局部最優(yōu)。3.實(shí)用性高。只需要一段描述日影變化過(guò)程的視頻以及相應(yīng)的時(shí)間就可以求出當(dāng)?shù)氐慕?jīng)緯度坐標(biāo)。模型的缺點(diǎn)：1.只能應(yīng)用在水平地面上。沒(méi)有考慮斜坡和坎坷地面上的情況。2.對(duì)于拍攝角度不好的視頻，該模型產(chǎn)生的誤差較大。3.對(duì)形狀規(guī)格較為復(fù)雜的建筑，本模型不能適用。九、 參考文獻(xiàn)[1][2]袁新生，邵大宏，LINGO和EXCEL在數(shù)學(xué)建模中的應(yīng)用，北京：科學(xué)出版社，2007。[3]謝金星，薛毅，優(yōu)化建模與LINDO/LINGO軟件，北京：清華大學(xué)出版社，[4]程序設(shè)計(jì)與應(yīng)用，高等教育出版社，2012。[5]孫榮峰，快速霍夫變換算法，計(jì)算機(jī)學(xué)報(bào)，24(10)：1104-1109，2001。[6]鄭鵬飛，林大鈞，基于影子軌線(xiàn)反求采光效果的技術(shù)研究，華東理工大學(xué)學(xué)報(bào)(自然科學(xué)版)，36(3)：458-463，2010。[7]何援軍，透視和透視投影變換，計(jì)算機(jī)輔助設(shè)計(jì)和圖形學(xué)學(xué)報(bào)，17(4)：734-739，2005。[8]武琳，基于太陽(yáng)陰影軌跡的經(jīng)緯度估計(jì)技術(shù)研究，天津大學(xué)，2010。[9]葉德意，基于自適應(yīng)變異粒子群算法的分布式電源選址與容量確定，電網(wǎng)技術(shù)，35(6)：155-157。十、 附錄10.1數(shù)據(jù)表格表11太陽(yáng)高度角與影長(zhǎng)數(shù)據(jù)時(shí)間太陽(yáng)高度角影長(zhǎng)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ù)時(shí)間節(jié)點(diǎn)影長(zhǎng)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源程序第一問(wèn)程序：clc;clear;shadow1=zeros(1,601);%初始化不同時(shí)間段影長(zhǎng)sthet=zeros(1,601);%初始化不同時(shí)間段太陽(yáng)高度角的sin值thet=zeros(1,601);%初始化不同時(shí)間段太陽(yáng)高度角n=295;%1月1日到10月22日間的天數(shù)time=[9:0.01:15];%從9點(diǎn)運(yùn)算到15點(diǎn)，步長(zhǎng)為0.01hh=[-48.60833333:0.15:41.39166667];%9點(diǎn)到15點(diǎn)地方恒星時(shí)系統(tǒng)下的時(shí)角delta=23.45*sin(2*pi*(284+n)/365);%太陽(yáng)赤緯height=3;%旗桿高度phi=39.9072222222222;%天安門(mén)緯度f(wàn)ori=1:601%從9點(diǎn)運(yùn)算到15點(diǎn)，共運(yùn)算601個(gè)值sthet(i)=cos(h(i)*pi/180)*cos(delta*pi/180)*cos(phi*pi/180)+sin(delta*pi/180)*sin(phi*pi/180);%計(jì)算出太陽(yáng)高度角的sin值thet(i)=asin(sthet(i))*180/pi;%計(jì)算出太陽(yáng)高度角shadow1(i)=sqrt((height/sthet(i))^2-(height)^2);%計(jì)算影長(zhǎng)endplot(time,shadow1)%輸出圖形ylabel('影長(zhǎng)/m');xlabel('時(shí)間/時(shí)');title('影子長(zhǎng)度與時(shí)間的關(guān)系');gridon;第二問(wèn)程序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; %北京時(shí)間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ù)氐胤綍r(shí)omega=15*(ts-12); %時(shí)角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));%太陽(yáng)高度角 a=2*pi-acos((sin(delta*pi/180)-sin(h)*sin(phi*pi/180))./(cos(h)*cos(phi*pi/180)));%太陽(yáng)方位角l=100./tan(h);%影長(zhǎng)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)%尋找最小點(diǎn)min=sum(abs(l-length)); %更新min為兩者更大的ganH=mean(Length.*tan(h));result=[lon(1,i)latganHabs(l-length)abs(a-angle)n];%最小點(diǎ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 %輸入容錯(cuò)error('error')endifnargout~=2 %輸出容錯(cuò)error('error')end%線(xiàn)形遞減策略w=MaxW-CurCount*((MaxW-MinW)/LoopCount);%得到粒子群群體大小以及一個(gè)粒子維數(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);%計(jì)算出全局最優(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;%計(jì)算每個(gè)粒子的新的適應(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ù))，進(jìn)行全局最優(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圖像處理程序%測(cè)試霍夫變換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);%計(jì)算二值圖像的標(biāo)準(zhǔn)霍夫變換，H為霍夫變換矩陣，theta,rho為計(jì)算霍夫變換的角度和半徑值figure,imshow(imadjust(mat2gray(H)),[],'XData',theta,'YData',rho,...'InitialMagnification','fit');xlabel('\theta(degrees)'),ylabel('\rho');axison,axisnormal,holdon;colormap(hot)%顯示霍夫變換矩陣中的極值點(diǎn)P=houghpeaks(H,50,'threshold',ceil(0.3*max(H(:))));%從霍夫變換矩陣H中提取5個(gè)極值點(diǎn)x=theta(P(:,2));y=rho(P(:,1));plot(x